{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load libraries\n",
    "import random\n",
    "import torch\n",
    "from plot_lib import set_default, show_mat\n",
    "from matplotlib.pyplot import plot, subplot\n",
    "from numpy import pi as π"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set style (needs to be in a new cell)\n",
    "set_default(figsize=(10, 4), dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate a collection of 5 3D vectors on a sphere\n",
    "A = torch.nn.functional.normalize(torch.randn(5, 3), p=2, dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking the content\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And the norm\n",
    "A.norm(p=2, dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change the value of the dimensionality d\n",
    "d = 8  # 1, 2, [3], 5, 8\n",
    "A = torch.nn.functional.normalize(torch.randn(10, d), 2, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this multiple times per given d\n",
    "a = A[[random.randrange(A.size(0))],:].t()  # perhaps comment, if playing with vis\n",
    "p = A @ a\n",
    "# p = torch.abs(A @ a)\n",
    "\n",
    "# Plot\n",
    "show_mat(A, a, p, threshold=-1)  # pick threshold ∈ [-1, +1), default -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the histogram of points uniformly distributed on a unit sphere in d dimensions\n",
    "bins = 101\n",
    "N = 10_000\n",
    "d_range = (2, 3, 4, 6, 10, 16); z = 1\n",
    "# d_range = d_range + (256,); z = 4  # uncomment\n",
    "ax1 = subplot(121)\n",
    "ax2 = subplot(122, projection='polar')\n",
    "for d in d_range:\n",
    "    A = torch.nn.functional.normalize(torch.randn(N, d), 2, 1)\n",
    "    B = torch.nn.functional.normalize(torch.randn(N, d), 2, 1)\n",
    "    h = torch.histc((A @ B.t()).view(-1), bins, -1, 1) / N**2\n",
    "    ax1.plot(torch.linspace(-1, 1, bins).numpy(), h.numpy())\n",
    "    ax2.plot(π * torch.linspace(-1, +1, 2 * bins).numpy(), torch.cat((h, h)).numpy() / 2)\n",
    "ax1.legend([f'{d}D' for d in d_range], ncol=3)\n",
    "ax1.set_ylim(0, 35e-3 * z)\n",
    "ax2.legend([f'{d}D' for d in d_range], ncol=3)\n",
    "ax2.set_rlim(0, 35e-3 * z / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These curves are simply the probability density function of a $\\Gamma$ distribution, with $\\alpha = \\beta = (D - 1) / 2$.\n",
    "See why this is the case [here](https://stats.stackexchange.com/a/85977/31844).\n",
    "More about curse of dimensionality and almost orthogonality [here](https://www.cs.princeton.edu/courses/archive/fall15/cos521/lecnotes/lec12.pdf)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:pDL] *",
   "language": "python",
   "name": "conda-env-pDL-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
